Kinetic energy cascades in quasi-geostrophic 
convection in a spherical shell 



Maxim Reshetnyak^, Pavel Hejda^ 

^ Institute of the Physics of the Earth, Russian Acad. Sci, 123995 Moscow, Russia 
E-mail: m . reshetnyakOgmail . com 

^ Institute of Geophysics, Academy of Sciences, 141 31 Prague, Czech RepubUc 
E-mail: ph@ig.cas.cz 

Abstract. We consider triadic nonlinear interaction in the Navier-Stokes equation for 
quasi-geostrophic convection in a spherical shell. This approach helps understanding 
the origin of kinetic energy transport in the system and the particular scheme of mode 
interaction, as well as the locality of the energy transfer. The peculiarity of convection 
in the sphere, concerned with excitation of Rossby waves, is considered. The obtained 
results are compared with our previous study in Cartesian geometry. 
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1. Introduction 

The redistribution of the physical fields is a quite widespread phenomenon in Nature. 
If we consider convective systems with dissipation, there are then two major effects 
which work in the opposite directions: diffusion (may be turbulent), which leads to 
the homogenization of the fields, and nonlinear interactions, which can produce large 
gradients and, in fact, cause the self-organization of the system. The most impressive 
examples of such behavior is the appearance of coherent structures in convective systems 
[T], topological pumping [2] and turbulent diamagnetism in the dynamo [3]. 

The considered effects demonstrate pumping of the fields in physical space to some 
domain of a scale smaller than the main scale of the system. If the sub domain is quite 
small, then it is clear that pumping should also be accompanied by the redistribution of 
the fields in wave space, e.g., the energy fiuxes between the different scales relevant to the 
nonlinear process, a-effect and separation in the scales of magnetic helicity [H |5] in the 
mean- field dynamo theory [6j, already mentioned coherent structures in convection. It is 
worth noting that fiuxes exist in the wave space even in quasistationary states provided 
that the scales of the energy injection into the system and the energy dissipation are 
different, see the introduction to the problem in [7] and the review of some recent results 
in [8j. 

In general, such fluxes can be deflned for various physical quantities, and the 
direction of their propagation is prescribed by many factors: by the nature of the 
considered quantity, the number of dimensions in the system and even the geometry 
of the system (as we shall see latter). Restricting our further research to convection, 
we recall that for the 3D homogeneous isotropic turbulence (in the absence of rotation) 
there is a direct energy cascade of the kinetic energy Ek from the large scales to the 
small scales, where dissipation takes place. The second invariant of the Navier-Stokes 
equation in 3D in the inviscid limit is kinetic helicity x- However, for isotropic turbulence 
its mean value is zero and its mean flux is also zero. 

The situation changes in 2D, where the cascade of kinetic energy is inverse: from 
the small to the large scales. Two-dimensional idealization was a useful approach for 
describing geophysical turbulence and was able to capture many important features of 
the flow. With regard to the more realistic models, one should consider quasi-geostrophic 
flows, which in view of their properties are somewhere inbetween 3D and 2D. In such a 
flow, one still has 3 dimensions, but the dependence of angular rotation {z coordinate) 
on direction is degenerated, see [9]. This flow is known by its structures elongated 
along z. The perpendicular scale of these structures is very small and deflned by the 
value of the Ekman number E. At the critical Rayleigh number, the flrst columnar 
mode has horizontal and vertical scales 0{E^^^), 0{1), respectively [10]. Note that 
E ~ 10~^^ in the Earth's liquid core, and the horizontal scale is extremely small. As 
follows from our previous study for the rotating rectangular box heated from below with 
periodic boundary conditions in the horizontal plane [11], these structures can supply 
kinetic energy in both directions: to the larger as well as to the smaller scales, which 
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can be important for understanding tlie energy budget in tlie planetary cores and for 
constructing semi-empirical models of turbulence. However, the solution of the same 
equations in plane geometry differs from that in the sphere. 

In Cartesian geometry, the increase of the heat sources characterized by the 
Rayleigh number leads to the gradual increase of the number of growing wave modes 
with smaller and larger scales than the scale of the first mode. However, in spherical 
geometry the increase of Ra excites Rossby waves. As a result the columns start to 
oscillate in the direction of the axis of rotation of the sphere, which leads to the rotation 
of columns around it. The direction of this rotation is defined by the slope of the outer 
boundary. For a spherical shell it is prograde and for a concave surface it has the 
opposite direction [12]. The appearance of this rotation (even differential) corresponds 
to the axi-symmetrical mode in spectral space. Thus the scenarios of energy transfer 
in the plane and in spherical geometries can be different. This is the motivation of 
our present study. Here, following ^1], we discuss the fluxes of kinetic energy in the 
standard Bousinesq equations used in the geodynamo for the different magnitudes of 
Ra and compare the results with the Cartesian geometry survey. 



2. Equations 

The thermal convection process driven by the flows of incompressible fluid (V ■ V = 0) 
in the Boussinesq approximation in a spherical shell (ricB < ^ < tcmb) rotating with 
the angular velocity Q in the z-direction is described by the Navier-Stokes equation 

Pr-^E(^^ + (V-V)V^ =-VP + F + eW (1) 

and the heat flux equation for temperature fluctuations T 
dT 

— + (V- V)(T + To) = V'T, (2) 

where To = cmb_ solution of the heat flux equation with flxed 

r CMB — ricB 

temperatures (1, 0) at the boundaries (rjcB, '''cmb ) in the absence of convection. 
Hereinafter ICB denotes the inner core boundary, CMB the core mantle boundary and 
(r, 9, if) is the spherical system of coordinates. 

The equations are scaled with the outer radius of the shell L, which makes the 
dimensionless radius tcmb = 1; the inner core radius tjcb is equal to 0.35, which is 
the value valid for the Earth. Velocity V, pressure P and the typical diffusion time t 
are measured in units of /t/L, p/t^/L^ and L^/k, respectively, where k is the thermal 
diffusivity, p is density, Pr = — is the Prandtl number, E = 2^\2 Ekman number 

and V is the kinematic viscosity. 

Force F includes the Coriolis and Archimedean effects: 



F = -1^ X V + RaTrl^, 



(3) 
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Figure 1. Meridional section of kinetic energy Ek f^or a) E = 2 10 ^, Pr — 1, Ra = 40, 
Rc - 17, (0, 60), with rotation, and b) E = 1, Pr = 1, Ra = 2.5 10^, Re ~ 3.4102, 
(0, 510^) without rotation (CorioHs force is switched off). Numbers in parentheses 
correspond to the range of the field. 



where Iz is the unit vector along the axis of rotation, Ra 



agoSTL 



is the modified 



Rayleigh number, a is the coefficient of volume expansion, 6T is the unit of temperature, 
and go is the gravitational acceleration at r = tcmb- 

The inner core, r < rjcB, with surface S'icb, can rotate around axis z due to the 
viscous torque, caused by the no-slip boundary conditions, used at both the boundaries 
'"iCB, ''"CMB- The momentum equation for the angular velocity u of the inner core has 
the form: 



where I is the inner-core moment of inertia along the z-axis. 

Equations (IT]-|1]) were solved using the standard spherical functions decomposition 
accompanied by a poloidal-toroidal decomposition of vector field V. The fast Chebyshev 
transform was used in the r-direction. The mesh grid in physical space was 128^. The 
Fortran code was parallelized in the r-direction using MPI. The details of the spherical 
function and Chebyshev polynomial decomposition can be found in [131 El HH] ■ 

3. Quasi-geostrophic convection 

Rapid rotation is a quite familiar phenomenon in geophysics. A strong Coriolis force 
leads to the appearance of elongated structures (columns) along the axis of rotation, 
in contrast to the cellular patterns in the non-rotating regime, where the spherically 
symmetrical (when the fields are averaged in time) buoyancy forces dominate, see Fig. [H 
These columns with a horizontal scale 1^ ~ E^^^ L ^ 1 rotate around their axes, so that 
the net helicity in the northern hemisphere for small Rayleigh numbers is negative and 
in the southern hemisphere positive, see Fig. |2]^a). Usually, the combination of the mean 




(4) 



SlCB 
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Figure 2. Meridional section of the mean kinetic helicity x for E = 2 10 ^, Pr = 1 a) 
Ra = 40, Re - 17, (-1.2 10^, 1.2 10^) and b) Ra = 4 10^, Re ~ 2 10^, (-3 10^ 3 10^). 
The dotted isohnes correspond to negative values. 

kinetic helicity x =< ^ " rot V > (closely connected with the a-effect) with differential 
rotation is considered as an explanation of the existence of the large-scale planetary 
magnetic fields. The helicity depends on the regularity of the flow and ratio of the 
Coriolis and Archemedean forces: the increase of the heat sources leads to the spread 
of convection to the Taylor cylinder (TC) accompanied by strong differential rotation. 
As a result, x changes sign in the middle of the spherical shell in TC, being positive at 
the inner-core boundary (ICB) and negative at the core-mantle boundary (CMB) in the 
northern hemisphere (the helicity in TC still retains its dipole structure with respect to 
the equatorial plane). Before we refer the reader to the more detailed study of kinetic 
helicity, including the influence of the magnetic field on x [iH], let us stress that, from 
the point of view of the mean- field theory , x in TC and outside it are not the same: 
in TC X is produced by the large-scale motions near the boundaries, and in the outer 
part it has a cyclonic nature with ^ 1, so that the separation of the fields into the 
large and small scales, adopted in the theory, is valid. 

4. Spectral properties 

As follows from the behavior of convection in physical space, rotation leads to the 
appearance of the small scales in the horizontal plane. Hereinafter we consider the 
spectral properties in the azimuthal direction, bearing in mind that the one- dimensional 
spectrum S{m) of field F{r, 6, (p) means: S{m) = J J f{m) f{m) r"^ sin 6 dr d6, where / 
is the Fourier transform of F and / is the complex conjugate of /. The spectra of three 
regimes are presented in Fig. [3l 
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m + 1 

Figure 3. Spectra of kinetic energy on the azimuthal wave number m with rotation for 
Ra — 40 (open circles), Ra = 800 (solid circles) and Ra = 2.5 10® (solid line) multiplied 
by 10^ without rotation. To use the logarithmic scale in the m-direction, the abscissa 
coordinate is shifted by 1. This means that the axi-symmetrical mode corresponds to 
value 1. 



Near the onset of convection there is the well-pronounced peak in the spectra 
(m = 3), corresponding to the horizontal scale of the column. Increase of Ra leads to 
the filling of the gaps in the spectra and transforming the saw-like structures to smooth 
curves. The slope of the spectrum for the developed rotating convection (solid circles) is 
less steep than for the non-rotating field because of the blocking of the nonlinear energy 
transfer over the spectrum by rotation. 

To study this problem we recall the technique discussed in [71 E] and our recent 
papers fU] [T7] . The main idea is to allocate a sphere of radius M in Fourier space 
and consider the evolution in time of kinetic energy E^{M) inside the sphere taking 
into account the integral flux of the kinetic energy n{M) from the outer part of the 
sphere in wave space with m > M. The exact form of the flux in physical space 
is n{M) = — [(V • V) V] ■ V^, where denotes the filtered velocity field with all 
harmonics with m > M equal to zero. As follows from [17], in Cartesian geometry 
the total flux n(A;), where k = |kx + ky + k^l, is mainly defined by the perpendicular 
contribution Il{k±) with k± = |kx + ky|. This is the motivation to consider further 
only the M-dependency. In contrast to our previous study [HI ITT] we do not proceed 

on 

with local- flux functions T = because the spherical case exhibits more irregular 

behavior of the flux than the flat geometry. Note also that U includes the local 
flux through a particular M, as well as the non-local contribution between the two 
subvolumes. 

All the curves, see Fig. H] for our regimes, satisfy condition 11 = for the large M 
conflrming the conservation of kinetic energy in the non-linear term of the Navier-Stokes 
equation. It is obvious that rotation changes the energy transfer substantially. In the 
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Figure 4. Integral flux Uk of kinetic energy as a function of the azimuthal 
wavenumber for the three regimes with the same labels as in FiglHl 

non-rotating regime, 11 is always negative, which means that the large scales feed the 
small scales (the direct cascade of energy). For the small M, the slope of increasing 11 is 
constant, which means that the flux to the small scales is non-local. The same situation 
holds for the larger M, where IJ increases, but the slope also remains constant. Due to 
this non-locality our solution differs from that for the Cartesian geometry in [HI [17] . 

For the geostrophic convection with small Ra we observe the inverse cascade of 
energy for M = — 2, n>0. The break of the curve at M = 3 means that, due to the 
sharp decrease of the spectrum, approximately the whole kinetic energy is in the sphere 
and n tends to zero. This is the reason why the direct cascade cannot be seen in the 
plot for large M. The increase of Ra decreases the relative flux to the axi-symmetrical 
mode twice. For larger M, 77 changes sign, which corresponds to the transition from 
the inverse cascade regime for M < 3 to the direct cascade of the kinetic energy for 
M > 3. 

The rotation changes not only the direction of energy transfer over the spectrum, 
but also the structure of the triangle in wave space, when two modes with m = P and 
m = Q produce the third mode m = K. To check this possibility, we have constructed 
the T2 flux function which describes the input of energy of the P-mode to the K-mode: 
T2{P, K) = - [(V(P) ■ V) V] ■ Y{K) Q. It is convenient to present T2 as a function of 
K — P, see Fig. |5l There is clear evidence of the non-local transfer of energy between 
the modes for small Ra under rotating convection. The shift of the extrema from the 
zero point describes the non-locality of interactions. This phenomenon can be explained 
as follows: the extrema correspond to the period (m = 3) between the local maxima in 
the spectra, see Fig. |3l Negative T2 for positive arguments means an inverse cascade. 
The increase of Ra leads to the complex interaction between the modes with both the 

J By definition T2 is antisymmetric function. 
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Figure 5. The fluxes of kinetic energy T2{K — P) for the three regimes, see Fig |3l4l 
Thin hne corresponds to the 5*'' order polynomial approximation. 

direct and inverse cascades with a different level of non-locality. The approximation of 
the curve with polynomials demonstrates the evident inverse cascade and non-locality 
in interactions, see Fig. [51 The comparison of the non-local characteristics for the 
spherical and Cartesian geometries demonstrates the increase of the non-locality transfer 
in the former case. There are many reasons for this phenomenon. Due to the spherical 
boundary, the spectra on m accumulate various scales in the 2;-directions. The pure 
statistics in the spherical problem may be more important: if the number of the columns 
is ~ k'j_ in the fiat geometry, their number is much smaller in the spherical geometry, 
because they are mainly distributed near TC. The other reason is that considered in 
[17] : k± = 8 for the columns is larger than that in the spherical geometry (M = 3) and 
we stand far away from the similarity region, even if we have similar resolutions in the 
models. That is why Cartesian geometry is so often used in MHD simulations to get 
scahngs. 

5. Discussion 

We have tried to present the regimes of geostrophic convection, well-known in the 
geodynamo community, by traditional tools of the turbulent community. Even if the 
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solution is quasi-stationary in physical space, there are non-zero fluxes of energy in wave 
space, proving that convection is not localized in wave space. Our results demonstrate 
that recent simulations for the usually used grids do not exhibit local transfer of energy, 
even without rotation. Rotation brings inverse flux into the first harmonics. In this way 
we can consider the excitation of the axi-symmetrical rotation of the columns around 
geographical axes to be the result of the inverse cascade. This phenomenon has no 
analog in the case of convection in the box, where rotation destroys rotating rolls with 
m = 0, for details of the Kiippers-Lortz instability refer to fi8\ [T9] . In some sense, the 
periodical boundary condition used in the Cartesian geometry corresponds to the case 
with the zero curvature of the boundary considered in [T2] . 

As can be seen from the presented simulations, the relative location of energy 
transfer will increase with increasing resolution of the model, but we believe that the 
mechanism of excitation of the axi-symmetrical flow due to the non-local inverse energy 
transfer will be the same. We hope that these results will also be interesting for solar 
modeling, where this technique could be used to check the energy exchange between 
the different radial layers and could thus help in constructing more realistic turbulent 
models. 
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